Export Genotypes for 9 Adaptive Run Timing Loci to NOAA ABL

Filtered for Plates Not Needing QC Crosschecks

Setup

rm(list = ls(all.names = TRUE))

if(!require("pacman")) install.packages("pacman"); library(pacman)

pacman::p_load(
  tidyverse,
  scales,
  janitor,
  GCLr,
  coin
)

knitr::opts_chunk$set(fig.width = 10)
knitr::opts_chunk$set(echo = TRUE)

.username = readLines("~/usr_pw.txt", n = 1)  # LOKI username
.password = readLines("~/usr_pw.txt" , n = 2)[[2]]  # LOKI password

Objective

The purpose of this notebook is to provide Wes Larson (NOAA ABL) with adaptive run timing genotypes (9 loci) for AHRP chum samples. Filter out plates undergoing QC crosschecks due to poor genotyping success. Wes is interested in representation across years, less across multiple streams. He wants to see how the patterns hold up over time. Might as well send everything we have though!

Steps

• import all existing genotypes for CM064-CM070
• filter out untrustworthy plates that need QC cross checks
• standard GT-seq genotype QA (missing, heterozygosity, duplicate)
• join sample date and spawning state <“../OceanAK/AHRP Salmon Biological Data 20240105_122404.368459_with_spawning_state.csv”>
• filter for 9 adaptive run timing loci
• export genotypes with date/spawning state

Background

An adaptive marker panel associated with run timing variation in chum salmon was developed by NOAA ABL and used to genotype two plates of samples from the ADFG Alaska Hatchery Research Program (AHRP) study being conducted within SE Alaska. AHRP samples (Fish Creek - Douglas Island 2023 + Prospect Creek 2017/2018) were genotyped in August 2024 for 20 single nucleotide polymorphisms and weighted linear regressions were performed to assess the association of different SNPs with run timing (collection date). The chum salmon adaptive marker panel contains 20 loci spread among 5 chromosomes and was developed from variation detected between Yukon River summer (Gisasa River) and fall (Pelly River) run chum salmon with low coverage whole genome sequencing. These populations spawn approximately two months and over 1500 river km apart. NOAA targeted the major peaks in differentiation with particular interest in markers within the genes leucine rich repeat containing nine (LRRC9; LG35) and estrogen receptor β (ESRB; LG29) as they were associated with return timing in other salmonid species. Eight SNPs (four on LG29 and four on LG35) appear to be associated with run timing within these collections. The strength of association for AHRP samples appears to vary by population. No markers on LGs 5, 23 or 15 appear to have a strong association with run timing within these collections. Overall, the results from LG29 and LG35 are promising and variation in these SNPs should be surveyed for a larger collection of samples.

A total of 9 adaptive, run timing loci were added to the AHRP SEAK Chum parentage panel:

• OkeV2_LG23_10546237
• OkeV2_LG29_25450752
• OkeV2_LG29_25455833
• OkeV2_LG29_25483595
• OkeV2_LG29_25515161
• OkeV2_LG35_28078867
• OkeV2_LG35_28128687
• OkeV2_LG35_28165076
• OkeV2_LG35_28167172

All AHRP Chum samples from Fish Creek - Douglas Island, Prospect Creek, Admiralty Creek, and Sawmill Creek were genotyped for the parentage panel at the ADF&G Gene Conservation Laboratory (GCL) in Winter/Spring 2025. Due to poor tissue quality, there was a high proportion of samples that failed to genotype and/or appear to be contaminated (high heterozygosity), resulting in low QC power for several plates. QC crosschecks are currently underway at the GCL in July 2025, however, Wes wants whatever genotypes we have on hand now so will just filter out those plates.

Import Data

LocusControl

Create objects, create LocusControl for Chum_AHRP_257_v1.3.0.

fs::dir_create("objects")
GCLr::create_locuscontrol(markersuite = "Chum_AHRP_257_v1.3.0",
                          username = .username,
                          password = .password)
GCLr::save_objects(objects = "LocusControl", path = "objects")

Hrmm, only 251 loci. Confirm that these are the same loci run in the most recent AHRP Chum project CM069

probe_loci_CM069 <- readr::read_tsv(
  file = "V:/Lab/Genotyping/SNP Projects/Chum/Project CM069 AHRP Parentage/Genotype Data Files/set2rr_CM069_outputs/probes.txt", 
  show_col_types = FALSE) %>%
  dplyr::distinct(`#Locus_ID`) %>%
  dplyr::pull()

all.equal(probe_loci_CM069, LocusControl$locusnames)

Yes, all are the same.

loci251 <- LocusControl$locusnames
GCLr::save_objects(objects = "loci251", path = "objects")

Genotypes

Read in genotypes by project to get sillyvec.

rm(LocusControl)
# GCLr::loki2r_proj(project_name = paste0("CM0", 64:70), username = .username, password = .password)  # comment out to avoid re-running

Confirm collections.

(sillyvec <- sort(project_sillys))

Subset sillyvec to remove alevin collection.

(sillyvec <- grep(pattern = "a", x = sillyvec, ignore.case = FALSE, value = TRUE, invert = TRUE))

Wipe out .gcl objects and re-read in using loki2r.

GCLr::save_objects(objects = "sillyvec", path = "objects")

rm(list = objects(pattern = ".gcl"))
rm(LocusControl)

Read in genotypes.

GCLr::load_objects(path = "objects", pattern = "LocusControl")

GCLr::loki2r(sillyvec = sillyvec,
             username = .username,
             password = .password, 
             test_type = "GTSNP")

Only took 5.39 minutes to read everything in at the office!

Save raw genotypes (useful for working offsite and off VPN).

fs::dir_create("data/genotypes")
fs::dir_create("data/genotypes/raw")
# GCLr::save_sillys(sillyvec = sillyvec, path = "data/genotypes/raw", rds = TRUE)  # hide so you don't overwrite on accident!

Re-load

Re-load, if necessary

GCLr::load_objects(path = "objects")

[1] "loci251"      "LocusControl" "sillyvec"    

GCLr::load_sillys(sillyvec = sillyvec, path = "data/genotypes/raw", rds = TRUE)

 [1] "CMADMCR13.gcl"   "CMADMCR14.gcl"   "CMADMCR17.gcl"   "CMADMCR18.gcl"   "CMFISHCR13.gcl"  "CMFISHCR14.gcl"  "CMFISHCR17.gcl"  "CMFISHCR18.gcl" 
 [9] "CMFISHCR19.gcl"  "CMFISHCR20.gcl"  "CMFISHCR21.gcl"  "CMFISHCR22.gcl"  "CMFISHCR23.gcl"  "CMFISHCRT13.gcl" "CMFISHCRT21.gcl" "CMFISHCRT22.gcl"
[17] "CMFISHCRT23.gcl" "CMPROSCR13.gcl"  "CMPROSCR14.gcl"  "CMPROSCR17.gcl"  "CMPROSCR18.gcl"  "CMPROSCR19.gcl"  "CMPROSCR20.gcl"  "CMPROSCR21.gcl" 
[25] "CMPROSCR22.gcl"  "CMPROSCRT21.gcl" "CMPROSCRT22.gcl" "CMSAWCR13.gcl"   "CMSAWCR14.gcl"   "CMSAWCR15.gcl"   "CMSAWCR17.gcl"   "CMSAWCR18.gcl"  
[33] "CMSAWCR19.gcl"   "CMSAWCR20.gcl"   "CMSAWCR21.gcl"   "CMSAWCR22.gcl"   "CMSAWCRT21.gcl"  "CMSAWCRT22.gcl" 

Metadata

Read in paired metadata from Salmon Biological Fact and Stream Specimens, format similar to pws_pink_ped_dat. NOTE I’m doing this in a bit of a hurry, so make sure to proof read carefully if recycling this code!

(metadata <- readr::read_csv(file = "../OceanAK/AHRP Salmon Biological Data 20240105_122404.368459_with_spawning_state.csv"))

Rows: 27046 Columns: 25
── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr  (10): silly_code, spawning_state, sex, tissue_type, dna_tray_code, location_code, sample_id, otolith_mark_present, otolith_mark_id, otolith_...
dbl   (9): collection_id, fish_id, length_mm, dna_tray_well_code, sample_year, is_missing_paired_data_exists, well_has_more_than_one_sample, fw_a...
lgl   (5): target_dna_tray_code, target_dna_tray_well_pos, target_container_array_type_id, container_array_type, determination_collection_id
date  (1): sample_date

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Modify to make similar to pws_pink_ped_dat, despite lack of location data from riverdist.

(
  metadata <- metadata %>% 
    dplyr::select(-dplyr::contains("target")) %>% 
    dplyr::select(-dplyr::contains("determination")) %>% 
    dplyr::rename(date = sample_date,
                  length = length_mm,
                  silly = silly_code) %>% 
    tidyr::unite(col = "sample", dna_tray_code:dna_tray_well_code, sep = "_", remove = FALSE) %>% 
    tidyr::unite(col = "silly_source", c(silly, fish_id), sep = "_", remove = FALSE) %>% 
    dplyr::mutate(franz_id = stringr::str_c(stringr::str_sub(silly, 1, 2),  # 1st letter of species name, 1st letter of stream name
                                            stringr::str_sub(silly, -2, -1),  # 2 digit sample year
                                            "_",
                                            stringr::str_pad(fish_id, width = 5, pad = "0", side = "left")),  # 5 digit fishid using padded zeros up front)
                  year = lubridate::year(date),
                  stream = stringr::str_remove(string = location_code, pattern = " Creek"),
                  DOY = lubridate::yday(date),  # day of year
                  sex = dplyr::case_when(
                    sex == "M" ~ "Male",
                    sex == "F" ~ "Female",
                    sex == "U" ~ NA_character_,
                    is.na(sex) ~ NA_character_),
                  origin = dplyr::case_when(
                    otolith_mark_present == "NO" ~ "Natural",
                    otolith_mark_present == "YES" ~ "Hatchery"
                  ),  # add origin variable
                  origin = base::factor(origin, levels = c("Natural", "Hatchery")),  # make factor to ensure hatchery != red
                  spawning_state = dplyr::case_when(
                    spawning_state == "Unknown" ~ NA_character_,
                    TRUE ~ spawning_state
                  )  # force Unknown to NA
    ) %>%  
    dplyr::select(
      franz_id,  # unique identifier for FRANz (e.g. PE16_00001, species, stream, year, GCL fish_id)
      sample,  # ultimate data key for AHRP, DWP barcode + DWP position
      silly_source,  # unique identifier for LOKI (GCL database)
      stream,  # 5 pedigree streams
      year,  # sample year (aka death year for FRANz)
      origin,  # hatchery or wild
      # hatchery,  # which hatchery
      sex,  # male, female, unknown
      date,  # sample date
      DOY,  # day of year (aka julian sample date)
      length,  # mideye to hypural plate length in mm
      # intertidal,  # sampling location above or below intertidal
      # distance_mouth,  # riverdist approximate sampling location from mouth in meters 
      spawning_state,  # alive, pink gill, grey gill, rotting (aka index of stream life)
      # pre_spawn,  # TRUE/FALSE (added in 2016); if the eggs and milt do not easily release and carcass has completely full gonads
      # partial_spawn,  # TRUE/FALSE (added in 2015); fish is sampled (dead or alive) with more than a little eggs or milt
      # preyed_upon,  # TRUE/FALSE (added in 2015); gonads removed by predators
      # stream_trib,  # stream tributary (Gilmour, Paddy, and Stockdale have multiple)
      # latitude,  # sampling area (where fish were moved to prior to sampling, rough approximate of spawning location)
      # longitude,  # sampling area (where fish were moved to prior to sampling, rough approximate of spawning location)
      otolith_mark_present,  # otolith mark present (yes = hatchery, no = natural)
      otolith_mark_id,  # otolith hatchery mark
      otolith_mark_status_code,  # otolith read status code (reason for no read = overground, lost, etc.)
      silly,  # GCL collection code
      fish_id,  # GCL fish ID within a silly
      dna_tray_code,  # 10-digit DWP barcode
      dna_tray_well_code,  # 48 DWP position 1:48
      # dna_tray_well_pos,  # 48 DWP position A1:H6
      # distance_tide,  # riverdist approximate sampling location from maximum high tide
      # riverdist_seg,  # riverdist river segment
      # riverdist_vert,  # riverdist vertex
      # riverdist_snapdist,  # riverdist snapping distance (aka how far was sampling location lat/long from the stream polyline)
      # high_tide  # highest extent of high tide in meters from mouth, unique to each stream
      fw_age,
      sw_age
    )  
)       

Sample Sizes

Verify that all genotyped samples have metadata.

geno_silly_source <- sapply(sillyvec, function(silly) {get(paste0(silly, ".gcl"))$SillySource}) %>% unlist()

table(geno_silly_source %in% metadata$silly_source)

FALSE  TRUE 
  241 18569 

Shoot, which samples do not have any metadata?

Ah, I bet they are the Fish Creek tag samples from 2013 that Tyler Dann and Heather Liller collected pre-AHRP field collections!

setdiff(geno_silly_source, metadata$silly_source)

  [1] "CMFISHCRT13_1"   "CMFISHCRT13_2"   "CMFISHCRT13_3"   "CMFISHCRT13_4"   "CMFISHCRT13_5"   "CMFISHCRT13_6"   "CMFISHCRT13_7"   "CMFISHCRT13_8"  
  [9] "CMFISHCRT13_9"   "CMFISHCRT13_10"  "CMFISHCRT13_11"  "CMFISHCRT13_12"  "CMFISHCRT13_13"  "CMFISHCRT13_14"  "CMFISHCRT13_15"  "CMFISHCRT13_16" 
 [17] "CMFISHCRT13_17"  "CMFISHCRT13_18"  "CMFISHCRT13_19"  "CMFISHCRT13_20"  "CMFISHCRT13_21"  "CMFISHCRT13_22"  "CMFISHCRT13_23"  "CMFISHCRT13_24" 
 [25] "CMFISHCRT13_25"  "CMFISHCRT13_26"  "CMFISHCRT13_27"  "CMFISHCRT13_28"  "CMFISHCRT13_29"  "CMFISHCRT13_30"  "CMFISHCRT13_31"  "CMFISHCRT13_32" 
 [33] "CMFISHCRT13_33"  "CMFISHCRT13_34"  "CMFISHCRT13_35"  "CMFISHCRT13_36"  "CMFISHCRT13_37"  "CMFISHCRT13_38"  "CMFISHCRT13_39"  "CMFISHCRT13_40" 
 [41] "CMFISHCRT13_41"  "CMFISHCRT13_42"  "CMFISHCRT13_43"  "CMFISHCRT13_44"  "CMFISHCRT13_45"  "CMFISHCRT13_46"  "CMFISHCRT13_47"  "CMFISHCRT13_48" 
 [49] "CMFISHCRT13_49"  "CMFISHCRT13_50"  "CMFISHCRT13_51"  "CMFISHCRT13_52"  "CMFISHCRT13_53"  "CMFISHCRT13_54"  "CMFISHCRT13_55"  "CMFISHCRT13_56" 
 [57] "CMFISHCRT13_57"  "CMFISHCRT13_58"  "CMFISHCRT13_59"  "CMFISHCRT13_60"  "CMFISHCRT13_61"  "CMFISHCRT13_62"  "CMFISHCRT13_63"  "CMFISHCRT13_64" 
 [65] "CMFISHCRT13_65"  "CMFISHCRT13_66"  "CMFISHCRT13_67"  "CMFISHCRT13_68"  "CMFISHCRT13_69"  "CMFISHCRT13_70"  "CMFISHCRT13_71"  "CMFISHCRT13_72" 
 [73] "CMFISHCRT13_73"  "CMFISHCRT13_74"  "CMFISHCRT13_75"  "CMFISHCRT13_76"  "CMFISHCRT13_77"  "CMFISHCRT13_78"  "CMFISHCRT13_79"  "CMFISHCRT13_80" 
 [81] "CMFISHCRT13_81"  "CMFISHCRT13_82"  "CMFISHCRT13_83"  "CMFISHCRT13_84"  "CMFISHCRT13_85"  "CMFISHCRT13_86"  "CMFISHCRT13_87"  "CMFISHCRT13_88" 
 [89] "CMFISHCRT13_89"  "CMFISHCRT13_90"  "CMFISHCRT13_91"  "CMFISHCRT13_92"  "CMFISHCRT13_93"  "CMFISHCRT13_94"  "CMFISHCRT13_95"  "CMFISHCRT13_96" 
 [97] "CMFISHCRT13_97"  "CMFISHCRT13_98"  "CMFISHCRT13_99"  "CMFISHCRT13_100" "CMFISHCRT13_101" "CMFISHCRT13_102" "CMFISHCRT13_103" "CMFISHCRT13_104"
[105] "CMFISHCRT13_105" "CMFISHCRT13_106" "CMFISHCRT13_107" "CMFISHCRT13_108" "CMFISHCRT13_109" "CMFISHCRT13_110" "CMFISHCRT13_111" "CMFISHCRT13_112"
[113] "CMFISHCRT13_113" "CMFISHCRT13_114" "CMFISHCRT13_115" "CMFISHCRT13_116" "CMFISHCRT13_117" "CMFISHCRT13_118" "CMFISHCRT13_119" "CMFISHCRT13_120"
[121] "CMFISHCRT13_121" "CMFISHCRT13_122" "CMFISHCRT13_123" "CMFISHCRT13_124" "CMFISHCRT13_125" "CMFISHCRT13_126" "CMFISHCRT13_127" "CMFISHCRT13_128"
[129] "CMFISHCRT13_129" "CMFISHCRT13_130" "CMFISHCRT13_131" "CMFISHCRT13_132" "CMFISHCRT13_133" "CMFISHCRT13_134" "CMFISHCRT13_135" "CMFISHCRT13_136"
[137] "CMFISHCRT13_137" "CMFISHCRT13_138" "CMFISHCRT13_139" "CMFISHCRT13_140" "CMFISHCRT13_141" "CMFISHCRT13_142" "CMFISHCRT13_143" "CMFISHCRT13_144"
[145] "CMFISHCRT13_145" "CMFISHCRT13_146" "CMFISHCRT13_147" "CMFISHCRT13_148" "CMFISHCRT13_149" "CMFISHCRT13_150" "CMFISHCRT13_151" "CMFISHCRT13_152"
[153] "CMFISHCRT13_153" "CMFISHCRT13_154" "CMFISHCRT13_155" "CMFISHCRT13_156" "CMFISHCRT13_157" "CMFISHCRT13_158" "CMFISHCRT13_159" "CMFISHCRT13_160"
[161] "CMFISHCRT13_161" "CMFISHCRT13_162" "CMFISHCRT13_163" "CMFISHCRT13_164" "CMFISHCRT13_165" "CMFISHCRT13_166" "CMFISHCRT13_167" "CMFISHCRT13_168"
[169] "CMFISHCRT13_169" "CMFISHCRT13_170" "CMFISHCRT13_171" "CMFISHCRT13_172" "CMFISHCRT13_173" "CMFISHCRT13_174" "CMFISHCRT13_175" "CMFISHCRT13_176"
[177] "CMFISHCRT13_177" "CMFISHCRT13_178" "CMFISHCRT13_179" "CMFISHCRT13_180" "CMFISHCRT13_181" "CMFISHCRT13_182" "CMFISHCRT13_183" "CMFISHCRT13_184"
[185] "CMFISHCRT13_185" "CMFISHCRT13_186" "CMFISHCRT13_187" "CMFISHCRT13_188" "CMFISHCRT13_189" "CMFISHCRT13_190" "CMFISHCRT13_191" "CMFISHCRT13_192"
[193] "CMFISHCRT13_193" "CMFISHCRT13_194" "CMFISHCRT13_195" "CMFISHCRT13_196" "CMFISHCRT13_197" "CMFISHCRT13_198" "CMFISHCRT13_199" "CMFISHCRT13_200"
[201] "CMFISHCRT13_201" "CMFISHCRT13_202" "CMFISHCRT13_203" "CMFISHCRT13_204" "CMFISHCRT13_205" "CMFISHCRT13_206" "CMFISHCRT13_207" "CMFISHCRT13_208"
[209] "CMFISHCRT13_209" "CMFISHCRT13_210" "CMFISHCRT13_211" "CMFISHCRT13_212" "CMFISHCRT13_213" "CMFISHCRT13_214" "CMFISHCRT13_215" "CMFISHCRT13_216"
[217] "CMFISHCRT13_217" "CMFISHCRT13_218" "CMFISHCRT13_219" "CMFISHCRT13_220" "CMFISHCRT13_221" "CMFISHCRT13_222" "CMFISHCRT13_223" "CMFISHCRT13_224"
[225] "CMFISHCRT13_225" "CMFISHCRT13_226" "CMFISHCRT13_227" "CMFISHCRT13_228" "CMFISHCRT13_229" "CMFISHCRT13_230" "CMFISHCRT13_231" "CMFISHCRT13_232"
[233] "CMFISHCRT13_233" "CMFISHCRT13_234" "CMFISHCRT13_235" "CMFISHCRT13_236" "CMFISHCRT13_237" "CMFISHCRT13_238" "CMFISHCRT13_239" "CMFISHCRT13_240"
[241] "CMFISHCRT13_241"

Correct. Add these to metadata with sample date info.

metadata <- 
  dplyr::bind_rows(metadata,
                   CMFISHCRT13.gcl %>% 
                     dplyr::select(SILLY_CODE, FK_FISH_ID, SillySource, DNA_TRAY_CODE, DNA_TRAY_WELL_CODE, CAPTURE_DATE) %>% 
                     dplyr::rename(silly = SILLY_CODE,
                                   fish_id = FK_FISH_ID,
                                   silly_source = SillySource, 
                                   dna_tray_code = DNA_TRAY_CODE,
                                   dna_tray_well_code = DNA_TRAY_WELL_CODE,
                                   date = CAPTURE_DATE) %>% 
                     tidyr::unite(col = "sample", dna_tray_code:dna_tray_well_code, sep = "_", remove = FALSE) %>% 
                     tidyr::unite(col = "silly_source", c(silly, fish_id), sep = "_", remove = FALSE) %>% 
                     dplyr::mutate(franz_id = stringr::str_c(stringr::str_sub(silly, 1, 2),  # 1st letter of species name, 1st letter of stream name
                                                             stringr::str_sub(silly, -2, -1),  # 2 digit sample year
                                                             "_",
                                                             stringr::str_pad(fish_id, width = 5, pad = "0", side = "left")),  # 5 digit fishid using padded zeros up front)
                                   date = lubridate::as_date(date),
                                   year = lubridate::year(date),
                                   stream = "Fish - Douglas Island",
                                   DOY = lubridate::yday(date)  # day of year
                     ) 
)

Verify all samples have metadata.

table(geno_silly_source %in% metadata$silly_source)

 TRUE 
18810 

Excellent, all present and accounted for.

Record which samples were genotyped.

metadata <- metadata %>% 
  dplyr::mutate(genotyped = silly_source %in% geno_silly_source)

How many samples genotyped per collection (silly)?

metadata %>% 
  dplyr::count(genotyped, silly) %>% 
  tidyr::pivot_wider(names_from = genotyped, values_from = n)

Check For Sample Duplicates

I got a weird many-to-many error in a join the first time I ran through this script, so I want to check and see if there are actually duplicate samples present and if so, remove them!

# bind all samples together into a single .gcl object
my.gcl <- sapply(sillyvec, function(silly) {
  my.silly <- get(paste0(silly, ".gcl"))
}, simplify = FALSE) %>% 
  dplyr::bind_rows()

Any duplicate SillySource?

table(table(my.gcl$SillySource))

    1 
18810 

Nope, all good there. Clean up.

Are there metadata duplicates?

table(table(metadata$silly_source))

    1     2 
27197    45 

Well son of a biscuit, there are!

metadata %>% 
  dplyr::filter(duplicated(silly_source) | duplicated(silly_source, fromLast = TRUE)) %>% 
  dplyr::arrange(silly_source)

Okay, they appear to be true duplicates in metadata, no conflicting columns. Remove duplicate rows.

metadata <- metadata %>% 
  dplyr::distinct(silly_source, .keep_all = TRUE)

Join Metadata

Join metadata to each silly by SillySource.

x <- sapply(sillyvec, function(silly) {
  my.silly <- get(paste0(silly, ".gcl"))
  my.silly <- my.silly %>% dplyr::left_join(
    y = metadata %>% dplyr::select(-genotyped),
    by = dplyr::join_by(SillySource == silly_source)
  )
  assign(x = paste0(silly, ".gcl"),
         value = my.silly,
         pos = 1)
})
rm(x)  # suppresses output

Filter QC Crosscheck Plates

Create a vector of plate IDs to toss, pending QC crosschecks.

plate_ids_qcxcheck <- c(64461, 64462, 64480, 64482, 64483, 64484, 64485, 64486, 64502, 64503, 64509, 64511, 64512, 64513, 64514, 64515, 64542, 64543, 64582, 56362, 56363) %>% as.character()

Get sample size by silly.

(
  sample_size_qcxcheck <- dplyr::tibble(silly = sillyvec) %>%
    dplyr::mutate(genotyped = GCLr::silly_n(sillyvec = sillyvec) %>% dplyr::pull(n))
)

Remove any fish on the plates pending QC crosschecks.

x <- sapply(sillyvec, function(silly) {
  my.silly <- get(paste0(silly, ".gcl"))
  my.silly <- my.silly %>% dplyr::filter(!stringr::str_detect(string = PLATE_ID, pattern = stringr::str_c(plate_ids_qcxcheck, collapse = "|")))
  assign(x = paste0(silly, ".gcl"),
         value = my.silly,
         pos = 1)
})
rm(x)  # suppresses output

How many fish were removed by silly?

(
  sample_size_qcxcheck <- sample_size_qcxcheck %>%
    dplyr::mutate(qcxcheck = genotyped - GCLr::silly_n(sillyvec = sillyvec) %>% dplyr::pull(n))
)

How many total samples?

sum(sample_size_qcxcheck$qcxcheck)

[1] 1995

Okay, that looks about right for 21 plates. Must have been some partial plates.

Genotype QA

Standard data QA (GT-seq, no conScore):

• Remove fish missing genotypes (<80% loci)
• Remove contamination (heterozygosity outliers +/- 3.5 modified z-score)
• Remove duplicates (>95% loci concordance, remove both duplicates)

(
  sample_size_qa <- dplyr::tibble(silly = sillyvec) %>%
    dplyr::mutate(genotyped = GCLr::silly_n(sillyvec = sillyvec) %>% dplyr::pull(n))
)

Missing Loci (<80%)

Remove fish with <80% loci genotyped.

miss_loci <-
  GCLr::remove_ind_miss_loci(sillyvec = sillyvec,
                             proportion = 0.8)

G3;30 IDs were removed from CMADMCR13.gcl
gG3;30 IDs were removed from CMADMCR17.gcl
gG3;4 IDs were removed from CMADMCR18.gcl
gG3;6 IDs were removed from CMFISHCR13.gcl
gG3;110 IDs were removed from CMFISHCR14.gcl
gG3;111 IDs were removed from CMFISHCR17.gcl
gG3;79 IDs were removed from CMFISHCR18.gcl
gG3;182 IDs were removed from CMFISHCR19.gcl
gG3;6 IDs were removed from CMFISHCR21.gcl
gG3;29 IDs were removed from CMFISHCR22.gcl
gG3;36 IDs were removed from CMFISHCR23.gcl
gG3;2 IDs were removed from CMFISHCRT13.gcl
gG3;2 IDs were removed from CMFISHCRT21.gcl
gG3;3 IDs were removed from CMFISHCRT22.gcl
gG3;1 IDs were removed from CMFISHCRT23.gcl
gG3;42 IDs were removed from CMPROSCR13.gcl
gG3;5 IDs were removed from CMPROSCR14.gcl
gG3;17 IDs were removed from CMPROSCR17.gcl
gG3;8 IDs were removed from CMPROSCR18.gcl
gG3;13 IDs were removed from CMPROSCR19.gcl
gG3;2 IDs were removed from CMPROSCR21.gcl
gG3;9 IDs were removed from CMPROSCR22.gcl
gG3;2 IDs were removed from CMPROSCRT22.gcl
gG3;23 IDs were removed from CMSAWCR13.gcl
gG3;1 IDs were removed from CMSAWCR14.gcl
gG3;6 IDs were removed from CMSAWCR15.gcl
gG3;46 IDs were removed from CMSAWCR17.gcl
gG3;17 IDs were removed from CMSAWCR18.gcl
gG3;17 IDs were removed from CMSAWCR19.gcl
gG3;8 IDs were removed from CMSAWCR21.gcl
gG3;25 IDs were removed from CMSAWCR22.gcl
gG3;A total of 872 individuals were removed from sillys in sillyvec.
g

# miss_loci$IDs_Removed

GCLr::save_objects(objects = "miss_loci", path =  "../objects", rds = TRUE)

(
  sample_size_qa <- sample_size_qa %>%
    dplyr::mutate(missing = genotyped - GCLr::silly_n(sillyvec = sillyvec) %>% dplyr::pull(n))
)

Overall 872 samples dropped due to poor quality DNA.

What was genotyping success by spawning state?

miss_loci_tib <- tibble::enframe(miss_loci$IDs_Removed, name = "silly", value = "fish_id") %>% 
  tidyr::unnest_longer(fish_id) %>% 
  dplyr::left_join(y = metadata, by = dplyr::join_by(silly, fish_id))

metadata %>% 
  dplyr::filter(genotyped) %>% 
  dplyr::mutate(drop_miss = silly_source %in% miss_loci_tib$silly_source) %>% 
  dplyr::count(drop_miss, spawning_state) %>% 
  tidyr::pivot_wider(names_from = drop_miss, values_from = n) %>% 
  dplyr::rename(n_drop = `TRUE`, n_keep = `FALSE`) %>% 
  dplyr::mutate(n = n_keep + n_drop,
                p_keep = n_keep / n,
                spawning_state = factor(spawning_state, levels = c("Alive", "Pink Gill", "Grey Gill", "Rotting"))) %>% 
  dplyr::arrange(spawning_state)

As expected, genotyping success was lower for rotting fish, but honestly, not terrible!

Contamination (Heterozygosity > +/- 3.5 modified Z-score)

Remove contaminated samples prior to checking for duplicates!

Individual heterozygosity outliers is our best metric for removing contaminated samples in the absence of something more sophisticated like conScore from GTscore. In previous analyses, including Shedd et al. 2022, we’ve used the standard 1.5 IQR outlier detection method to remove excessively heterozygous individuals. However, the 1.5 IQR method assumes a normal distribution, but heart sample heterozygosities tend to be skewed right due to contamination. As of 2024-02-08, we decided to pivot to using modified Z-scores with cutoffs of +/- 3.5 as recommended by Iglewicz and Hoaglin ¹. The modified Z-score uses the median and median absolute deviation (MAD) instead of the mean and standard deviation, and is thus more robust to outliers and asymmetrical distributions.

Source function and calculate individual heterozygosities from AHRP-PWS-Pink-Salmon-Pedigrees GitHub repo.

source("~/GitHub_repos/AHRP-PWS-Pink-Salmon-Pedigrees/code/functions/calc_ind_het.R")

(
  ind_het <-
    calc_ind_het(
      sillyvec = sillyvec,
      loci = loci251,
      ncores = parallel::detectCores() - 4
    )
)

Time difference of 17.24236 secs

Calculate +/- 3.5 modified z-score cutoffs by stream. NOTE this means that the exact heterozygosity cutoff will vary slightly by stream. For PWS Pink Salmon, I did it separately by silly, but some of our sample sizes are pretty small here and we have the separate “tag” sillys. I’m opting to do this by stream (all years together) to allow for subtle changes in allele frequency per stream.

# Function to calculate median absolute deviation (MAD)
mad <- function(x) {
  median(abs(x - median(x)))
}

# Calculate the modified z-score for each individual
(
  ind_het <- ind_het %>%
    dplyr::mutate(
      year = 2000 + as.numeric(stringr::str_sub(
        string = silly,
        start = -2,
        end = -1
      )),
      stream = dplyr::case_when(
        stringr::str_detect(string = silly, pattern = "ADM") ~ "Admiralty",
        stringr::str_detect(string = silly, pattern = "FISH") ~ "Fish",
        stringr::str_detect(string = silly, pattern = "PROS") ~ "Prospect",
        stringr::str_detect(string = silly, pattern = "SAW") ~ "Sawmill",
        TRUE ~ "mistakes_were_made"
      )
    ) %>%
    dplyr::group_by(stream) %>%  # cutoffs are specific to each stream!
    dplyr::mutate(
      modified_z_score = 0.6745 * (het - median(het)) / mad(het),
      outlier = dplyr::case_when(abs(modified_z_score) > 3.5 ~ TRUE, TRUE ~ FALSE)
    ) %>%
    dplyr::ungroup()
)

  
GCLr::save_objects("ind_het", path = "objects", rds = TRUE)

Plot distribution of heterozygosity and show outliers.

ind_het %>% 
  ggplot2::ggplot(ggplot2::aes(x = het, fill = outlier)) +
  ggplot2::geom_histogram(binwidth = 1 / length(loci251)) +
  ggplot2::facet_grid(rows = ggplot2::vars(stream), scales = "free_y") +
  ggplot2::xlim(0, 1) +
  ggplot2::scale_fill_manual(name = "Outlier", values = c("TRUE" = "black", "FALSE" = "grey60")) +
  ggplot2::xlab("Individual Heterozygosity") +
  ggplot2::ylab("Frequency") +
  ggplot2::ggtitle("Individual Heterozygosity - By Stream") +
  ggplot2::theme_bw(base_size = 14) 

Remove outliers.

output <- lapply(sillyvec, function(x) {
  GCLr::remove_ids(
    silly = x,
    IDs = ind_het %>% dplyr::filter(outlier, silly == x) %>% dplyr::pull(fish_id)
  )
})

G3;8 IDs were removed from CMADMCR13.gcl
gG3;5 IDs were removed from CMADMCR14.gcl
gG3;147 IDs were removed from CMADMCR17.gcl
gG3;2 IDs were removed from CMADMCR18.gcl
gG3;11 IDs were removed from CMFISHCR13.gcl
gG3;113 IDs were removed from CMFISHCR14.gcl
gG3;111 IDs were removed from CMFISHCR17.gcl
gG3;49 IDs were removed from CMFISHCR18.gcl
gG3;55 IDs were removed from CMFISHCR19.gcl
gG3;2 IDs were removed from CMFISHCR20.gcl
gG3;7 IDs were removed from CMFISHCR21.gcl
gG3;17 IDs were removed from CMFISHCR22.gcl
gG3;100 IDs were removed from CMFISHCR23.gcl
gG3;1 IDs were removed from CMFISHCRT13.gcl
gG3;0 IDs were removed from CMFISHCRT21.gcl
gG3;1 IDs were removed from CMFISHCRT22.gcl
gG3;3 IDs were removed from CMFISHCRT23.gcl
gG3;18 IDs were removed from CMPROSCR13.gcl
gG3;10 IDs were removed from CMPROSCR14.gcl
gG3;153 IDs were removed from CMPROSCR17.gcl
gG3;6 IDs were removed from CMPROSCR18.gcl
gG3;50 IDs were removed from CMPROSCR19.gcl
gG3;0 IDs were removed from CMPROSCR20.gcl
gG3;2 IDs were removed from CMPROSCR21.gcl
gG3;10 IDs were removed from CMPROSCR22.gcl
gG3;0 IDs were removed from CMPROSCRT21.gcl
gG3;0 IDs were removed from CMPROSCRT22.gcl
gG3;17 IDs were removed from CMSAWCR13.gcl
gG3;2 IDs were removed from CMSAWCR14.gcl
gG3;4 IDs were removed from CMSAWCR15.gcl
gG3;93 IDs were removed from CMSAWCR17.gcl
gG3;1 IDs were removed from CMSAWCR18.gcl
gG3;0 IDs were removed from CMSAWCR19.gcl
gG3;0 IDs were removed from CMSAWCR20.gcl
gG3;2 IDs were removed from CMSAWCR21.gcl
gG3;0 IDs were removed from CMSAWCR22.gcl
gG3;0 IDs were removed from CMSAWCRT21.gcl
gG3;0 IDs were removed from CMSAWCRT22.gcl
g

message(paste0("A total of ", ind_het %>% dplyr::filter(outlier) %>% nrow(), " individuals were removed from sillys in sillyvec."))

G3;A total of 1000 individuals were removed from sillys in sillyvec.
g

Not too many contaminated samples! Most likely culled during missing loci.

(
  sample_size_qa <- sample_size_qa %>%
    dplyr::mutate(heterozygosity = genotyped - missing - GCLr::silly_n(sillyvec = sillyvec) %>% dplyr::pull(n))
)

Overall 1,000 samples dropped due to sample contamination.

What was contamination by spawning state?

metadata %>% 
  dplyr::filter(genotyped) %>% 
  dplyr::mutate(het_outlier = silly_source %in% (ind_het %>% dplyr::filter(outlier) %>% dplyr::pull(SillySource))) %>% 
  dplyr::count(het_outlier, spawning_state) %>% 
  tidyr::pivot_wider(names_from = het_outlier, values_from = n) %>% 
  dplyr::rename(n_drop = `TRUE`, n_keep = `FALSE`) %>% 
  dplyr::mutate(n = n_keep + n_drop,
                p_keep = n_keep / n,
                spawning_state = factor(spawning_state, levels = c("Alive", "Pink Gill", "Grey Gill", "Rotting"))) %>% 
  dplyr::arrange(spawning_state)

As expected, contamination (heterozygosity outliers) was higher for rotting fish, but honestly, not terrible!

Duplicate (>=95%)

Identify potential duplicate genotypes (>=95% loci).

duplicate_check_95 <-
  GCLr::dupcheck_within_silly(
    sillyvec = sillyvec,
    minproportion = 0.95,
    minnonmissing = 0.6,
    ncores = parallel::detectCores() - 4
  )

Time difference of 21.62875 secs

GCLr::save_objects("duplicate_check_95", path = "objects", rds = TRUE)

How many duplicate pairs by silly?

duplicate_check_95 %>% 
  dplyr::count(silly) %>% 
  dplyr::arrange(dplyr::desc(n))

Remove both duplicates! As opposed to GSI work, where we want to keep individuals but aren’t typically worried about paired data, here we want to remove both individuals as the paired data integrity (including sample date, otolith reads, etc.) is lost.

duplicates_removed <-
  GCLr::remove_dups(dupcheck = duplicate_check_95, remove_both = TRUE)

G3;2 IDs were removed from CMADMCR14.gcl
gG3;20 IDs were removed from CMADMCR17.gcl
gG3;6 IDs were removed from CMADMCR18.gcl
gG3;2 IDs were removed from CMFISHCR13.gcl
gG3;12 IDs were removed from CMFISHCR14.gcl
gG3;14 IDs were removed from CMFISHCR17.gcl
gG3;38 IDs were removed from CMFISHCR18.gcl
gG3;18 IDs were removed from CMFISHCR19.gcl
gG3;4 IDs were removed from CMFISHCR23.gcl
gG3;2 IDs were removed from CMFISHCRT23.gcl
gG3;2 IDs were removed from CMPROSCR13.gcl
gG3;4 IDs were removed from CMPROSCR14.gcl
gG3;14 IDs were removed from CMPROSCR17.gcl
gG3;4 IDs were removed from CMPROSCR18.gcl
gG3;10 IDs were removed from CMPROSCR19.gcl
gG3;4 IDs were removed from CMPROSCR20.gcl
gG3;2 IDs were removed from CMPROSCR22.gcl
gG3;8 IDs were removed from CMSAWCR17.gcl
gG3;4 IDs were removed from CMSAWCR18.gcl
gG3;2 IDs were removed from CMSAWCR22.gcl
g

GCLr::save_objects("duplicates_removed", path = "objects", rds = TRUE)

(
  sample_size_qa <- sample_size_qa %>%
    dplyr::mutate(
      duplicate = genotyped - missing - heterozygosity - GCLr::silly_n(sillyvec = sillyvec) %>% dplyr::pull(n)
    )
)

Overall 157 samples dropped due to duplicate genotypes.

Final

Final, post-QA sample sizes by silly.

(
  sample_size_qa <- sample_size_qa %>%
    dplyr::mutate(
      final =  GCLr::silly_n(sillyvec = sillyvec) %>% dplyr::pull(n)
    )
)

GCLr::save_objects("sample_size_qa", path = "objects", rds = TRUE)

fs::dir_create("output")
readr::write_csv(x = sample_size_qa, file = "output/sample_size_qa.csv")

Save

Save final genotypes and update metadata.

Genotypes

Save post-QA genotypes with joined metadata.

fs::dir_create("data/genotypes/postQA")
GCLr::save_sillys(sillyvec = sillyvec, path = "data/genotypes/postQA", rds = TRUE)

Metadata

add pass_QA

postQA_silly_source <- sapply(sillyvec, function(silly) {get(paste0(silly, ".gcl"))$SillySource}) %>% unlist()

metadata <- metadata %>% 
  dplyr::mutate(pass_QA = silly_source %in% postQA_silly_source)

Confirm.

metadata %>% 
  dplyr::count(genotyped, pass_QA)

NOTE that genotyped in sample_size_qa doesn’t match genotyped in metadata because I pre-filtered out QC Crosscheck plates.

sample_size_qa %>% 
  janitor::adorn_totals() %>% 
  dplyr::filter(silly == "Total")

 silly genotyped missing heterozygosity duplicate final
 Total     16815     872           1000       157 14786

QA Conclusions

Overall 18,810 total pedigree samples were genotyped; 1,995 were removed due to pending QC Crosschecks; 872 were removed due to poor quality DNA (missing >= 20% loci); 1,000 were removed due to contaminated DNA (heterozygosity >= 3.5 modified z-score); and 157 were removed due to duplicate genotypes (>= 95% loci), leaving a total of 14,786 in the final dataset.

Adaptive Loci Genotypes

Create a single object with just the metadata and adaptive loci genotypes for NOAA ABL.

Calculate Allele Frequencies

I’m going to look at this three different ways:

1. naive, all samples (including hatchery strays) with collection date,
2. natural-origin fish only (no hatchery strays) with collection date, and
3. natural-origin fish only (no hatchery strays) with adjusted collection date based on spawning_state (re-scaling to Alive) * Alive = collection date * Pink Gill = collection date - 2 days * Grey Gill = collection date - 5 days * Rotting = collection date - 10 days

Stream sampling was limited by both crew availability and weather conditions (flooding), so collection date is not a perfect proxy for run timing. Additionally, in some of the later years crews were tagging live fish near the intertidal to boost genetic sample size; so these live fish were sampled much earlier than post-spawn carcasses with the same run timing.

All Samples (Hatchery + Natural)

Calculate the allele counts for each locus by stream, year, and DOY. Adapted from Pat Barry’s code in AHRP_ChumAdaptive.html

loci_adaptive <- grep(pattern = "OkeV2", x = loci251, value = TRUE)
GCLr::save_objects(objects = "loci_adaptive", path = "objects/")

# testing
# chum_adaptive_geno %>%
#   dplyr::select(silly_source, stream, year, DOY, tidyselect::starts_with(loci_adaptive[1])) %>%
#   tidyr::pivot_longer(cols = tidyselect::starts_with(loci_adaptive[1]), names_to = "locus", values_to = "allele") %>%
#   dplyr::mutate(allele = factor(allele, levels = LocusControl$alleles[[loci_adaptive[1]]]$call),
#                 allele = as.numeric(allele)) %>%  # make factor and convert to numeric
#   dplyr::count(stream, year, DOY, allele) %>%
#   dplyr::mutate(locus = loci_adaptive[1])

# calculate the number of each allele observed per stream, year, day
# NOTE: there is no filtering for hatchery/natural origin
(chum_adaptive_allele_counts <- lapply(loci_adaptive, function(locus) {
  chum_adaptive_geno %>%
    dplyr::select(silly_source,
                  stream,
                  year,
                  DOY,
                  tidyselect::starts_with(locus)) %>%
    tidyr::pivot_longer(
      cols = tidyselect::starts_with(locus),
      names_to = "locus",
      values_to = "allele"
    ) %>%
    dplyr::mutate(allele = factor(allele, levels = LocusControl$alleles[[{{locus}}]]$call),  # order alleles as per LocusControl
                  allele = as.numeric(allele)) %>%  # make factor and convert to numeric
    dplyr::count(stream, year, DOY, allele) %>%
    dplyr::mutate(locus = locus)
}) %>%
    dplyr::bind_rows())

Convert to allele frequencies.

(
  chum_adaptive_allele_freq <- chum_adaptive_allele_counts %>%
    tidyr::pivot_wider(names_from = allele, values_from = n) %>%
    dplyr::select(-`NA`) %>%  # drop no-calls
    dplyr::mutate(
      `1` = replace_na(`1`, 0),
      # account for unobserved alleles on a given day
      `2` = replace_na(`2`, 0)
    ) %>%
    tidyr::pivot_longer(
      cols = `1`:`2`,
      names_to = "allele",
      values_to = "n"
    ) %>%
    dplyr::group_by(stream, year, DOY, locus) %>%
    dplyr::mutate(
      n_alleles = sum(n),
      p_alleles = n / n_alleles,
      SNP = paste(locus, allele, sep = "_")
    ) %>% 
    dplyr::ungroup() %>% 
    dplyr::filter(n_alleles > 0)  # remove any stream/year/DOY/locus combinations without data
)

Confirm data structure.

chum_adaptive_allele_freq %>% 
  dplyr::group_by(locus) %>% 
  dplyr::summarise(min(p_alleles),
                   max(p_alleles))

chum_adaptive_allele_freq %>% 
  dplyr::filter(is.na(p_alleles))

Natural-origin Samples (no Hatchery Strays)

Re-calculate the allele counts for natural-origin fish only for each locus by stream, year, and DOY. No hatchery strays! Adapted from Pat Barry’s code in AHRP_ChumAdaptive.html

# calculate the number of each allele observed per stream, year, day
# NOTE: removing all hatchery-origin strays!!!
(chum_adaptive_allele_counts_wild <- lapply(loci_adaptive, function(locus) {
  chum_adaptive_geno %>%
    dplyr::filter(origin == "Natural") %>%  # no hatchery strays or unknown origin fish
    dplyr::select(silly_source,
                  stream,
                  year,
                  DOY,
                  tidyselect::starts_with(locus)) %>%
    tidyr::pivot_longer(
      cols = tidyselect::starts_with(locus),
      names_to = "locus",
      values_to = "allele"
    ) %>%
    dplyr::mutate(allele = factor(allele, levels = LocusControl$alleles[[{{locus}}]]$call),  # order alleles as per LocusControl
                  allele = as.numeric(allele)) %>%  # make factor and convert to numeric
    dplyr::count(stream, year, DOY, allele) %>%
    dplyr::mutate(locus = locus)
}) %>%
    dplyr::bind_rows())

Convert to allele frequencies.

(
  chum_adaptive_allele_freq_wild <- chum_adaptive_allele_counts_wild %>%
    tidyr::pivot_wider(names_from = allele, values_from = n) %>%
    dplyr::select(-`NA`) %>%  # drop no-calls
    dplyr::mutate(
      `1` = replace_na(`1`, 0),
      # account for unobserved alleles on a given day
      `2` = replace_na(`2`, 0)
    ) %>%
    tidyr::pivot_longer(
      cols = `1`:`2`,
      names_to = "allele",
      values_to = "n"
    ) %>%
    dplyr::group_by(stream, year, DOY, locus) %>%
    dplyr::mutate(
      n_alleles = sum(n),
      p_alleles = n / n_alleles,
      SNP = paste(locus, allele, sep = "_")
    ) %>% 
    dplyr::ungroup() %>% 
    dplyr::filter(n_alleles > 0)  # remove any stream/year/DOY/locus combinations without data
)

Adjust for Spawning State

Re-calculate the allele counts for natural-origin fish only for each locus by stream, year, and DOY. No hatchery strays! Adjust DOY based on spawning state, re-scaling to Alive. Adjustments are verrrry rough guesses.

• Rotting = Alive + 10 days
• Grey Gill = Alive + 5 days
• Pink Gill = Alive + 2 days

Adapted from Pat Barry’s code in AHRP_ChumAdaptive.html

# calculate the number of each allele observed per stream, year, day
# NOTE: removing all hatchery-origin strays!!!
(chum_adaptive_allele_counts_wild_adj_spawning <- lapply(loci_adaptive, function(locus) {
  chum_adaptive_geno %>%
    dplyr::filter(origin == "Natural") %>%  # no hatchery strays or unknown origin fish
    dplyr::mutate(DOY = dplyr::case_when(spawning_state == "Pink Gill" ~ DOY - 2,
                                         spawning_state == "Grey Gill" ~ DOY - 5,
                                         spawning_state == "Rotting" ~ DOY - 10,
                                         TRUE ~ DOY
                                         )) %>%  # adjust DOY based on spawning date, adjustments are complete WAGs!!!
    dplyr::select(silly_source,
                  stream,
                  year,
                  DOY,
                  tidyselect::starts_with(locus)) %>%
    tidyr::pivot_longer(
      cols = tidyselect::starts_with(locus),
      names_to = "locus",
      values_to = "allele"
    ) %>%
    dplyr::mutate(allele = factor(allele, levels = LocusControl$alleles[[{{locus}}]]$call),  # order alleles as per LocusControl
                  allele = as.numeric(allele)) %>%  # make factor and convert to numeric
    dplyr::count(stream, year, DOY, allele) %>%
    dplyr::mutate(locus = locus)
}) %>%
    dplyr::bind_rows())

Convert to allele frequencies.

(
  chum_adaptive_allele_freq_wild_adj_spawning <- chum_adaptive_allele_counts_wild_adj_spawning %>%
    tidyr::pivot_wider(names_from = allele, values_from = n) %>%
    dplyr::select(-`NA`) %>%  # drop no-calls
    dplyr::mutate(
      `1` = replace_na(`1`, 0),
      # account for unobserved alleles on a given day
      `2` = replace_na(`2`, 0)
    ) %>%
    tidyr::pivot_longer(
      cols = `1`:`2`,
      names_to = "allele",
      values_to = "n"
    ) %>%
    dplyr::group_by(stream, year, DOY, locus) %>%
    dplyr::mutate(
      n_alleles = sum(n),
      p_alleles = n / n_alleles,
      SNP = paste(locus, allele, sep = "_")
    ) %>% 
    dplyr::ungroup() %>% 
    dplyr::filter(n_alleles > 0)  # remove any stream/year/DOY/locus combinations without data
)

View Allele Frequencies

Create functions to standardize plots. Each set of plots is facetted to show years/loci. NOTE since I am forcing the y-axes to 0-1, the SE ribbons get cut off for some of the fit lines.

adaptive_freq_plot <- function(freq, stream) {
  freq %>% 
    dplyr::filter(stream == {{stream}}) %>% 
    dplyr::mutate(locus = stringr::str_remove_all(string = locus, pattern = "OkeV2_")) %>%  # shorten locus names to fit
    ggplot2::ggplot(ggplot2::aes(x = DOY, y = p_alleles, colour = SNP, group_by(SNP))) +
    ggplot2::geom_point() +
    ggplot2::geom_smooth(method = "lm", weight = "n_alleles", formula = y~x) +
    ggplot2::ylim(0, 1) +  # this cuts off the SE's for the fit line
    ggplot2::facet_grid(locus ~ year) + 
    ggplot2::theme_bw() + 
    ggplot2::theme(legend.position = 'none',
                   axis.text.x = ggplot2::element_text(angle = 90, vjust = 0.5, hjust = 0.5)) +
    ggplot2::labs(title = stream, x = "DOY", y = "Allele Frequency")
}

avg_allele_n_plot <- function(freq, stream) {
  freq %>% 
    dplyr::filter(stream == {{stream}}) %>% 
    dplyr::group_by(stream, year, DOY) %>% 
    dplyr::summarise(n_alleles = round(mean(n_alleles)), .groups = "drop") %>% 
    ggplot2::ggplot(ggplot2::aes(x = DOY, y = n_alleles)) + 
    ggplot2::geom_col() +
    ggplot2::facet_grid(~year) +
    ggplot2::labs(title = paste0(stream, " - Average Allele Sample Size per Day"), y = "Average Number of Alleles") +
    ggplot2::theme_bw() +
    ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90, vjust = 0.5, hjust = 0.5))
}

All Samples

Includes hatchery-origin strays and unknown origin fish (i.e. no otolith read).

Admiralty

Show average allele sample sizes per DOY by year.

avg_allele_n_plot(freq = chum_adaptive_allele_freq, stream = "Admiralty")

NOTE small sample sizes in 2018 and later sampling dates due to flooding conditions.

adaptive_freq_plot(freq = chum_adaptive_allele_freq, stream = "Admiralty")

Inconsistent patterns among years for most loci. The most consistent relationship is in LG29_25515161. Some evidence of a brood year effect with more similarities between 2013/2017 (assuming most chum are age 4).

Fish

Show average allele sample sizes per DOY by year.

avg_allele_n_plot(freq = chum_adaptive_allele_freq, stream = "Fish - Douglas Island")

adaptive_freq_plot(freq = chum_adaptive_allele_freq %>% dplyr::filter(year != 2020), stream = "Fish - Douglas Island")  # ditch 2020

Inconsistent patterns among years for LG35 loci, but LG29 loci are reasonably consistent.

Prospect

Show average allele sample sizes per DOY by year.

avg_allele_n_plot(freq = chum_adaptive_allele_freq, stream = "Prospect")

NOTE sample sizes are lackluster for 2020 and 2021, as are the tails of the run in 2018.

adaptive_freq_plot(freq = chum_adaptive_allele_freq, stream = "Prospect")

Inconsistent patterns among years for LG29 loci (see LG29_25455833 for 2014 vs. 2019!), but LG35 loci are reasonably consistent.

Sawmill

Show average allele sample sizes per DOY by year.

avg_allele_n_plot(freq = chum_adaptive_allele_freq, stream = "Sawmill")

NOTE sample sizes reasonable for all except 2020.

adaptive_freq_plot(freq = chum_adaptive_allele_freq %>% dplyr::filter(year != 2020), stream = "Sawmill")  # ditch 2020; low n

More similar to Fish Creek - Douglas Island. Weak, inconsistent relationship for LG35 loci, but LG29 loci have a consistent relationship across years.

Natural-origin

Only includes natural-origin fish.

Admiralty

Show average allele sample sizes per DOY by year.

avg_allele_n_plot(freq = chum_adaptive_allele_freq_wild, stream = "Admiralty")

Very few hatchery strays into Admiralty, no big changes in sample sizes.

adaptive_freq_plot(freq = chum_adaptive_allele_freq_wild, stream = "Admiralty")

Similar patterns as noted above.

Fish

Show average allele sample sizes per DOY by year.

avg_allele_n_plot(freq = chum_adaptive_allele_freq_wild %>% dplyr::filter(year != 2020), stream = "Fish - Douglas Island") # ditch 2020; low n

Lot of hatchery strays in Fish Creek - Douglas Island! Good sample sizes throughout the run for most years, excluding 2013 (sparse sampling) and 2021 (low n).

adaptive_freq_plot(freq = chum_adaptive_allele_freq_wild %>% dplyr::filter(year != 2020), stream = "Fish - Douglas Island")  # ditch 2020; low n

Some subtle changes in relationships when excluding hatchery strays, but general pattern of stronger relationship with LG29 loci holds.

Prospect

Show average allele sample sizes per DOY by year.

avg_allele_n_plot(freq = chum_adaptive_allele_freq_wild, stream = "Prospect")

NOTE sample sizes remain lackluster for 2020 and 2021, as are the tails of the run in 2018.

adaptive_freq_plot(freq = chum_adaptive_allele_freq_wild, stream = "Prospect")

Minimal changes to relationships after excluding hatchery strays.

Sawmill

Show average allele sample sizes per DOY by year.

avg_allele_n_plot(freq = chum_adaptive_allele_freq_wild, stream = "Sawmill")

NOTE sample sizes reasonable for all except 2020.

adaptive_freq_plot(freq = chum_adaptive_allele_freq_wild %>% dplyr::filter(year != 2020), stream = "Sawmill")  # ditch 2020

The pattern of stronger relationships with LG29 loci holds.

Natural-origin + Adjust for Spawning State

Only includes natural-origin fish and adjustments to collection date (DOY) based on spawning state (re-scaled to Alive).

Admiralty

Show average allele sample sizes per DOY by year.

avg_allele_n_plot(freq = chum_adaptive_allele_freq_wild_adj_spawning, stream = "Admiralty")

Distribution changes a bit due to spawning_state adjustment to DOY.

adaptive_freq_plot(freq = chum_adaptive_allele_freq_wild_adj_spawning, stream = "Admiralty")

Adjusting DOY for spawning_state doesn’t appear to alter relationships much for Admiralty. LG29 loci have a stronger association with run timing than LG35.

Fish

Show average allele sample sizes per DOY by year.

avg_allele_n_plot(freq = chum_adaptive_allele_freq_wild_adj_spawning %>% dplyr::filter(year != 2020), stream = "Fish - Douglas Island") # ditch 2020; low n

Lot of hatchery strays in Fish Creek - Douglas Island! Good sample sizes throughout the run for most years, excluding 2013 (sparse sampling) and 2021 (low n).

adaptive_freq_plot(freq = chum_adaptive_allele_freq_wild_adj_spawning %>% dplyr::filter(year != 2020), stream = "Fish - Douglas Island")  # ditch 2020

Adjusting DOY for spawning_state definitely introduces some more noise, but LG29 loci still have a stronger association with run timing than LG35.

Prospect

Show average allele sample sizes per DOY by year.

avg_allele_n_plot(freq = chum_adaptive_allele_freq_wild_adj_spawning, stream = "Prospect")

NOTE bigger shift in the distribution of DOY for 2013 and 2019 than other years.

adaptive_freq_plot(freq = chum_adaptive_allele_freq_wild_adj_spawning, stream = "Prospect")

Adjusting DOY for spawning_state definitely introduces some more noise, but cleans up the relationships with LG29 loci (except for LG29_25455833, dunno what is going on with those brood year effects!).

Sawmill

Show average allele sample sizes per DOY by year.

avg_allele_n_plot(freq = chum_adaptive_allele_freq_wild_adj_spawning, stream = "Sawmill")

NOTE sample sizes getting fairly small per DOY, but good distribution.

adaptive_freq_plot(freq = chum_adaptive_allele_freq_wild_adj_spawning %>% dplyr::filter(year != 2020), stream = "Sawmill")  # ditch 2020

Adjusting DOY for spawning_state appears to have reversed the relationships for 2019 for LG29_2540752 and LG_25483595?

Export Genotypes for NOAA ABL

Save as both a .csv and .rds.

readr::write_csv(x = chum_adaptive_geno, file = "output/chum_adaptive_geno_2025-07-30.csv")
GCLr::save_objects(objects = "chum_adaptive_geno", path = "output", rds = TRUE)

Conclusions

Nine adaptive loci associated with run timing were identified by NOAA ABL and included in the ADF&G AHRP chum salmon parentage GT-seq panel. ~19K SEAK chum AHRP pedigree samples were genotyped at the GCL. ~4K samples were removed from further analyses due to pending QC cross checks, poor quality DNA, contaminated samples, or duplicate genotypes, leaving a total of ~15K samples collected from 2013-2023. Weighted linear regressions were performed to investigate the association between adaptive loci allele frequencies and run timing (collection date) using code adapted from Pat Barry (NOAA ABL). Note that collection date is an imperfect proxy of run timing for these fish given that samples were collected from both living fish and post-spawn carcasses! Adjusting collection date for fish condition (spawning_state) appears to have improved some relationships. The strength of association appears to vary both by stream, year, and origin (hatchery vs. natural), with some apparent brood year effects (assuming most chum are age-4). The sole LG23 SNP shows inconsistent and at times contrasting patterns across years, whereas SNPs on LG29 (ESRB) and LG35 (LRRC9) show more consistency across streams and years. LG29 loci (except LG_25455833) were more consistently associated with run timing than LG35 loci in all streams after adjusting DOY to account for spawning_state.

end

---

1. Boris Iglewicz and David Hoaglin (1993), “Volume 16: How to Detect and Handle Outliers”, The ASQC Basic References in Quality Control: Statistical Techniques, Edward F. Mykytka, Ph.D., Editor.↩︎